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Abstract 



Oh! We estimate the QCD coupling constant from a lattice calculation of the 



bottomonium spectrum. The second order perturbative expansion of the pla- 



k> \ quette expectation value is employed to determine as at a scale set by the 

H : 

>y , 2S-1S and IP-IS level splittings. The latter are computed in NRQCD in 

a dynamical gauge field background with two degenerate flavours of Wilson 
quarks at intermediate masses and extrapolated to the chiral limit. Com- 
bining the Nf = 2 result with the quenched result at equal lattice spacing 
we extrapolate to the physical number of light flavours to find a value of 
a^jJmz) = 0.1118(17). The error quoted covers both statistical and system- 
atic uncertainties in the scale determination. An additional 5% uncertainty 
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comes from the choice of the underlying sea quark formulation and from trun- 
cation errors in perturbative expansions. 

I. INTRODUCTION 

A precise knowledge of the strong coupling constant is of central importance both for 
strong-interaction phenomenology and for stringent tests of the standard model as a whole. 
There has been significant progress in determining as and its running during the last years: 
Measurements are available for a large number of different reactions covering energies up to 
190 GeV. Although 'the' global average has remained nearly unchanged the uncertainties of 
single measurements as well as the scatter between them has been strongly reduced. The 
clustering of results into groups of high-energy and low-energy data which was observed 
back in 1995 has completely disappeared and the various results now nicely align around 
an average of C(y^{mz) ~ 0.119 |[l|]. While the majority of determinations of the strong 
coupling relies on perturbative expansions of high-energy cross sections, lattice simulations 
of QCD provide, in principle, a nonperturbative access to determine as from low-energy 
quantities (a review of lattice methodology is presented in 0). In fact, a viable scheme 
along these lines has been demonstrated with the Schrodinger functional technique in the 
context of the quenched approximation, Nf = |^. An alternative approach is to consider 
perturbative expansions of suitable short- distance quantities and to determine their vacuum 
expectation values from simulations. In this manner Monte Carlo data can be utilized to 
extract estimates for as at a scale which is provided by the lattice spacing, a 0]. This 
procedure lends itself easily to the setting of full QCD, given that sufficient samples of 
QCD vacuum configurations are available as demonstrated some time ago in ref. 0. In the 
meantime, the statistical precision of the resulting lattice estimates of as in Nf = and 2 
theories has been increased such as to allow for an extrapolation in the number of dynamical 
flavours to Nf = 3. 

Although, a priori, any mass or energy splitting calculated in a lattice simulation can be 



used to set the lattice spacing, a, in practice one clearly prefers those choices which are least 
prone to systematic errors. Spin-averaged radial and orbital angular momentum splittings 
in heavy quarkonium bound states, in particular in bottomonium, are favourite candidates 
. They are largely insensitive to the heavy quark mass: experimentally one finds the spin- 
averaged mass splittings between IP and IS levels as well as 2S and IS levels to be nearly 
equal for the T and \1' bound states. Moreover, one may infer phenomenologically, that 
they depend only mildly on the light quark masses: assuming characteristic gluon momenta 
exchanged between quarks inside the T to be O(lGeV), heavy flavours will be negligible in 
virtual quark loops whereas the masses of light quarks are equally small compared to this 
momentum scale so that their exact values as well as the differences between non-strange 
and strange quark masses are likely to be of little importance. These arguments imply 
that (a) there is no need to tune the bare heavy quark mass precisely and (b) to a good 
approximation the contribution of light quarks can be mimicked by considering three light 
flavours of average mass ^. The latter feature is a major beneflt in lattice simulations 
which are technically restricted to rather large quark masses and generally rely on chiral 
extrapolations. 

The bonus in choosing spin-averaged level splittings in bottomonium to set the scale 
stands against the malus that the direct simulation of heavy quark dynamics employing 
standard relativistic actions suffers from serious discretization errors. Among several alter- 
natives to circumvent this problem, the nonrelativistic effective theory for QCD (NRQCD) 
PI provides the most appropriate technique for simulating systems containing a b-quark. 
The NRQCD Lagrangian is written as a series of operators ordered according to powers of 
the mean-squared heavy-quark velocity. Each operator introduces a new coupling constant 
which is determined by perturbatively matching the theory to QCD. The decoupling of 
quark and antiquark flelds in NRQCD is equivalent to integrating out the heavy quark mass 
and, as a consequence, the amount of lattice spacing errors is governed by the characteris- 
tic momentum exchanged between heavy quarks, ©(aAqco); rather than by the underlying 
quark mass scale. 



The realistic inclusion of vacuum polarization effects constitutes a serious hurdle in to- 
day's lattice simulations of QCD, in particular for Wilson fermions. One should remember, 
that the state-of-the-art updating algorithms require an even number of degenerate sea 
quarks. As a consequence the strong coupling cannot be determined on field configurations 
with three kinds of active sea quarks, u, d, s; rather one has to resort to an extrapolation in 
the number of active flavours. Using the staggered discretisation for the dynamical fermions, 
Davies and coworkers p presented a careful analysis of 0:5 (including the various sources 
of error). Consistent results, although with substantially larger errors, have been quoted by 
the authors of [|],|lOl who also use staggered sea quarks but heavy Wilson valence quarks 
instead of nonrelativistic quarks. 

In this paper, we shall apply the methods of ||^ to the case of dynamical Wilson fermions. 
They carry different flnite-a errors (compared to staggered fermions) that affect both the 
nonperturbative results for the bottomonium splittings and the plaquette expectation value 
as well as the perturbative expansion of the plaquette. We shall improve on previous error 
analyses by studying the quark mass dependence of quarkonium splittings; hence we shall 
be able to reduce the uncertainty due to the light quark mass scale. 

The paper is organized as follows. In Sec. |I| we briefly recall the set-up of SESAM's 
hybrid Monte Carlo simulation of lattice QCD with dynamical Wilson fermions as well as 
the relevant issues of the computation of nonrelativistic propagators from a lattice NRQCD 



action. We proceed with a presentation of spectroscopy results in Sec. |T|. Since within 
the T system the singlet states rjf, and hb have not been observed in experiment yet, we 
have to deviate from the suggestion to use fully spin-averaged splittings and have to choose 
2^5*1 — 1'^ 5*1 and ^P — l^Si to determine the lattice spacing. Here, ^P denotes the spin average 
of ^Pj. Therefore the discussion of flt results and systematic errors in Sec. ^ covers radial 
splittings as well as spin splittings. In Sec. ^ we determine the quenched and unquenched 
plaquette couplings and perform the extrapolation in Nf. We conclude in Sec. |V| with a 
discussion of the resulting value of as in the continuum MS scheme. 



II. SIMULATION SET-UP 



The present work is based on ref. [11 and can be described as the final analysis of 
SESAM lattices with increased statistics, both with respect to the number of configurations 
and NRQCD propagator inversions, to reduce the error. 

A. Gauge Fields 

We have performed a hybrid Monte Carlo simulation of full QCD at /3 = 5.6 with 
two degenerate fiavours of dynamical standard Wilson fermions. A single lattice size of 
L^ X T = 16^ X 32 is used. It corresponds to a physical box of 1.2-1.4 fm in the spatial 
direction (depending on the experimental quantity used to set the scale) which is sufficiently 
large to exclude finite volume effects on the bottomonium ground state and the first radial 
and orbital angular momentum excitation. We generated configurations at three different 
values of the sea quark hopping parameter, n, each sample consisting of 5000 trajectories 
from which 200 decorrelated vacuum configurations are chosen. Table | lists the simulation 
parameters. For details of the HMC run and subtleties concerning autocorrelation studies we 
refer the reader to [|T^ . The quenched run uses an overrelaxed Cabbibo-Marinari heat bath 
update, thermalised with 2000 sweeps, and measurements are performed on configurations 
separated by 250 sweeps. 

B. Some NRQCD Prerequisites 



Compared to our intermediate results presented in ||I3|, we have significantly increased 
our statistics at k = 0.1575 and k = 0.1570 corresponding to the two lighter quark masses. 
Heavy quark propagators are calculated in the nonrelativistic approximation. We have im- 
plemented the NRQCD action at 'next-to-leading order' with spin-independent operators 
of 0(r?7,6f^), 0{mhV^) and spin-dependent terms of 0(r7ibf^) and Olm^v^) included. The 
radiative corrections to the coefficients of these interaction terms are not exactly accounted 



for. Rather we use the mean field prescription to cover their main effect and choose the 
tadpole parameter computed from the mean link in Landau gauge. Unlike for light hadrons, 
in lattice simulations of heavy quarkonium bound states one can efficiently increase statis- 
tics by computing propagators from different, widely separated source points on the same 
configuration. Hence we have averaged the data from up to 12 point sources for the Nf = 2 
lattices. Note, that we do not use a multi-source, but evolve each starting point separately. 
In the quenched runs we restrict ourselves to 4 source points as new configurations can be 
generated with comparatively little effort. To project on radially excited and orbital angular 
momentum states some care is needed in choosing an approriate smearing procedure. The 
method of choice for bottomonium is a fixed-gauge wavefunction smearing technique. For 
the T and rif, we calculate a 4 x 4 matrix of correlators with four different smearings at 
source and sink, sc/sk = 1, 1, 2, 3, corresponding to a point source (1), the ground state (1), 
the first (2) and second (3) excited states, respectively. For the L = 1 states we restrict 
ourselves to the ground state and the first excitation as signals deteriorate. The smearing 



functions that we use are the wavefunctions calculated in Ref. [0|. With all the coefficients 
in the action set to their tree level values the only parameter left besides the gauge cou- 
pling is the bare heavy quark mass, ami,. We did not tune the bare b-quark mass but kept 
arub = 1.7 throughout the simulation. This value reproduces the correct T mass in the 
quenched approximation ||13| and it turns out to be adequate in the full theory, too, leading 
to kinetic masses mki„(T) = 9.97(28), 9.63(24), 9.68(27) GeV for k = 0.1560,0.1570,0.1575, 
respectively. 

III. BOTTOMONIUM SPECTROSCOPY 

A. Fit Results 

We determine the triplet-S state l^Si and its first radial excitation as well as the singlet- 
P ground state l^Pi by a simultaneous fit of two source-smeared correlators to a double- 



exponential ansatz 

Csc,iT) = bl,e-^^^^ + bi,e-'^^^^ sc= 1,2. (1) 

These fits yield the cleanest signals and are very stable. We varied the fit interval over 
a considerable range as illustrated in Fig. ^ for the sample with the lightest sea quark 
content. The analysis reveals clear plateaus in global masses for each meson correlator. 
Since fits to smeared-local correlators usually tend to overestimate the plateau value, we 
have varied the fitting procedure to check for the stability of the results. First, we have 
used correlators which are smeared both at source and sink in a two-exponential fit. Sink 
smearing produces noisier signals. Nevertheless, ground state energies are obtained with 
similar accuracy and yield consistent results. Excited states, however, exhibit significantly 
larger errors compared to the smeared-local data. Second, we have applied fits involving 
two smeared-local correlators but three exponentials. The third exponential is meant to 
account for contaminations from higher radial excitations and to provide a more reliable 
estimate of the 2S state. The results are in very good agreement with those obtained from 
two exponentials, thus providing confidence in the radially excited level. Since both singlet 
resonances, r]f, and hh, have not been observed experimentally we have to recourse to the 
known splittings 2'^S'i — 1^5*1 and P — l^Si to set the lattice scale. Here, P denotes the spin- 
average triplet-P state: ^P = | {^Pq + 3 ■'^ Pi + 5 ■'^ P2). To obtain ^P we have to determine 
the P fine structure which is accomplished by single exponential fits to the ratio of two 
correlators 

C,{T)/C2{T) = A exp{-AEuT) , AEu = ^1 - P2 • (2) 

This way we compute the splittings of ^Pj relative to ^Pi for J = 0, 1,2 which are then 
combined to form the spin average. In the case of smeared-local correlators the plateau 
sets in only at rather large times. For P-wave states these can hardly be reached before 
the signal is drowned into noise. Results quoted for spin splittings have therefore been 
taken from smeared-smeared correlators. These run into plateaus at early times where they 



differ considerably from their smeared-local counterparts as illustrated in Fig. |T|. Table |I| 
summarizes the relevant fit results. Errors are taken from 300 bootstraps. 

B. Extrapolation in the light quark mass 

The hopping parameters that we have chosen correspond to quark masses in the range 
ms/2 < niq < rris where rris denotes the strange quark mass. According to ||15| one expects 



the energy splittings to depend linearly on the light quark mass, mg, in lowest order of a 
chiral expansion. Hence we choose the simple ansatz 

aAE = aAEo + c anig (3) 

to extrapolate in rrig. As illustrated in Fig. |^ the splittings considered here exhibit a very 
mild dependence on the sea quark mass. In fact, one cannot distinguish their values at 
am = arris/S ~ 0.0159 and our smallest simulated mass.[| As outlined in the introduction m 
is a reasonable average quark mass provided the strange-nonstrange mass difference can be 
neglected for the lowest T bound states. If extrapolated down to the physical light quark 
mass the splittings increase modestly. The parameters of uncorrelated linear fits are listed in 
Table |T| together with the splittings at rh and mi. We conclude that within the accuracy of 
our data no deviation from linear behaviour is found and we emphasize that sea quark mass 
dependences are smaller by one order of magnitude compared to the light hadron sector |]12 



We may choose the difference in mean values of aAE{m,s/3) and aAE{mi) as an upper 
limit of the uncertainty in the sea quark mass: l^P — 1^5*1 exhibits a 4% decrease, whereas 
2^Si — I'^S'i is affected by 7%. From the experimental bottomonium splittings T' — T = 
0.5629 GeV and x — T = 0.4398 GeV one determines the lattice scales summarized in Table 
[rV| . We use the average value of a~^ at rh to convert our lattice results into physical units. 



^Note that Kg determined using the K mass disagrees with the value calculated from K* or $. 
We use an average here since the exact value of Kg is clearly not important in the present analysis. 
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As intended, it matches the (mean) value of the quenched lattice spacing at /3 = 6.0, so that 
unquenching effects can be studied. 

It is obvious from Fig. ^that the quenched spectrum does not reproduce the experimental 
spectrum while the Nf = 2 data is in much closer agreement. This behaviour can be ascribed 
to the different running of the coupling in both theories and translates into a discrepancy 
between the plaquette coupling as determined from the 2S-1S and IP- IS splittings in the 
quenched approximation as we shall see below. 

C. Error Estimates 

Besides the dependence on the sea quark mass several error sources deteriorate the 
accuracy of the energy splittings and hence affect the precision in scale determination: 

Inherent to the nonrelativistic approach one is faced with higher order relativistic cor- 
rections due to the truncation of the effective NRQCD Lagrangian of choice and to the 
incomplete renormalization of coefficients in the velocity expansion. Having included spin- 
independent terms through 0{mf)V^) we expect that spin-averaged splittings will receive a 
1% correction from higher orders (interactions of 0{mf,v^)). But this is just the order of 
magnitude incurring by second-order spin-dependent effects! Hence we may estimate the 
truncation error by switching off the latter. By inspection of Table we find that this 
procedure leaves the energy levels unaltered. A ratio fit to the S hyperfine splitting, how- 
ever, reveals a 15% decrease when 0{mbV^) interactions are switched on as is expected from 
power counting. The impact of the incomplete renormalization of the NRQCD expansion 
coefficients is exposed by changing the tadpole improvement factor from the Landau-link 



prescription, Uq = (|Trf/^)LG, to the plaquette prescription, Uq = ^ {^TrU^„). It turns 
out to be of similar magnitude. As a conservative estimate, we shall henceforth quote an 
additional error of 0.10 in a^^ to cover both effects and refer to it as the uncertainty in the 
NRQCD expansion. 

We have explicitly calculated the dependence of the 2S-1S radial splitting on the bare 



heavy quark mass. Of course, we expect very little change as we vary amt. This is indeed 
confirmed as can be seen from Fig. ^, where we plot the triplet S radial splitting as a 
function of kinetic mass. The latter is computed for each bare mass value, arrib = 1.6 — 2.0, 
by giving the meson a small amount of momentum and fitting the nonrelativistic dispersion 
relation. 

We have not checked for finite volume effects. We expect them to be much smaller for 
quarkonia than for light mesons. Quenched lattices of size ~ 1.5 fm are found to be sufficient 
for the lowest charmonium levels [01. Hence our results for 2S-1S and IP-IS in bottomonium 



are safe. But, as pointed out in |T^, higher radial excitations like 2P or 3S require a linear 
lattice extent larger than 2 fm, even for bottomonium. 

Finally, our results will be affected by various finite-a errors. Note, that these cannot 
be removed by extrapolation, since there exists no continuum limit to NRQCD. Instead 
the latter should be viewed as being based on an effective action that is geared to obtain 
physical results at finite cut-off only. Obviously, in oder to make sure that the effective 
action approach is a useful tool one has to ascertain that the spectrum is independent of 
the lattice spacing within a certain window. This then would establish that the matching 
of the NRQCD action to QCD is sufficiently accurate. 

For nonrelativistic b-quarks discretization errors are likely to be larger than in the light 
hadron spectrum. The improvement of the NRQCD Lagrangian is thus a crucial issue and 
its efficiency has to be checked explicitly by simulating at different values of f3. We have 
removed Ola"^) errors in the tree approximation, but we cannot perform a scaling analysis 
for the dynamical data, as we are restricted to a single lattice spacing. In the quenched 



approximation Davies et. al. [|T7| have studied the bottomonium spectrum on three quenched 
lattices with spacings in the range 0.05 fm to 0.15 fm. They find good scaling in the ratio 
of radial and orbital bb splittings to the p-mass if the latter is computed from a tadpole 
improved clover action. Also ratios of such splittings within the T-system do not exhibit 
any dependence on the lattice spacings (spin splittings do!). Although these results are 
quite encouraging, we emphasize that dynamical Wilson fermions introduce additional linear 
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scaling violations whose size can only be safely estimated through simulations at different 
lattice spacings. 

IV. PLAQUETTE COUPLING 

A value of the strong coupling now is readily obtained: compute the expectation value 
of a short distance quantity on the lattice and match it with the perturbative expansion. 
Obviously, this is a reasonable procedure only, if nonperturbative effects are negligible which 
is definitely the case for the simplest lattice quantity, the 1x1 Wilson loop. While the 
common choice of the expansion parameter in the continuum is the MS-coupling, on the 
lattice it is more suitable to choose a subtraction scheme that refers to a nonperturbative 
quantity like the static QQ-potential |18[. Rather than using ay itself we adopt here a 



slightly modified scheme, defined in Ref. through 

- \Yi{hlTiU^,) = yap {^^ [1 - (1.1870 + 0.0249 Nj) ap] , (4) 

the rationale behind this definition of plaquette coupling, ap, being a matter of convenience: 
one has to worry about higher order perturbative corrections only once, when converting 
the lattice coupling into a standard continuum scheme at the very end of the analysis. Note 
that Eq. | is valid in the chiral limit and for Wilson fermions. One prefers to expand the 
logarithm of the plaquette since it converges more rapidly than the plaquette expectation 
value itself. The scale 3.41/a is the 'average gluon momentum' in the first-order contribution 



to — Iniyii computed with the technique suggested in ||T9[. In Table we summarise the 
couplings ap obtained from Eq.(^ as well as the scales determined from the IP-IS {x " S) 
and 2S-1S (T' - T) splittings. 

In the unquenched case we quote values for both arriq = arris/ 3 and arriq = ami to 
estimate the systematic error connected to the finite sea-quark mass. Plaquette expectation 
values have been extrapolated accordingly. We do not quote an error for them since it is 
negligible compared to the uncertainty in the scale. Subsequently these couplings are evolved 
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to a common scale /i = 9.0 GeV using the universal two-loop P function, see Table [VII| . The 
error in the evolution is minute as the evolution range is very small. For instance, using just 
the one-loop evolution results in a deviation of much less than 1%. 

The plaquette couplings in the quenched and unquenched theories can now be extrap- 
olated to the number of active light quark flavours which is expected to be Nj = 3 in the 
case of the low-lying bb bound states. Guided by the perturbative evolution, we extrapolate 
ap^ linearly in Nj, Figure ^. Obviously the mismatch between ap- values obtained from dif- 
ferent splittings in the quenched approximation disappears, once the dynamical quarks are 
switched on. Thus we find, that full QCD with Wilson fermions shows similar behaviour as 
QCD with staggered dynamical fermions 0]. Nevertheless, there is no way to pin down the 
the number of active flavours precisely. Continued extrapolation to Nf = 4 further increases 
ap (with the error blown up) and leads to la increase in a^^lMz). Hence simulations with 
Nf > 2 are called for to reduce this uncertainty. 

Consider once more the ordering of extrapolations that have been performed to arrive 

(3) 

at ap . First we have extrapolated splittings in lattice units as a function of the dynamical 
quark mass to determine a~^ at m. Then we performed the A^/-extrapolation. One might 
argue about this prescription of carrying out the limits. Alternatively one might prefer to 
do the flavour extrapolation at fixed quark mass. Although in principle at finite lattice 
spacing the order may matter, we see no difference: the couplings separately obtained at 



each sea-quark mass, see Table VIII, do not reveal any significant dependence on m„ and 



are consistent with those in Table VII 



V. DISCUSSION 

To make the connection with the MS-scheme one invokes 

»m (Q) = ^p'' i^'^'Q) [i + l^p'' + C2mi<y'p''? + 0((«?^^¥ 

with a scale factor e~^/^ chosen to eliminate the Nf dependence in the first-order coefficient 
of the expansion pO). The crucial point about Eq.(^ concerns the coefficient C2 which is 
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(5) 



only known in pure gauge theory and thus causes a significant uncertainty on a|^. Let us, for 
the moment, ignore this uncertainty and set C2 ~ 0.95, as in the quenched theory, to study 
the errors directly related to the lattice method. We start from Op (9.0 GeV) and evolve the 
coupling to the Z meson mass scale, applying the formulae in |^: First, a^(e~^/^-9.0GeV), 
is evolved down to the charm threshold with the three-loop beta function. Matching the 
three flavour with the four flavour theory, one obtains a^{Mc) which in turn is evolved 
upwards to the b-quark threshold. We have chosen ttIc = 1.3 GeV and mj, = 4.1 GeV for 



the charm and bottom thresholds, see Table |^. As was already noted in p2| , the value of 
the coupling at irtz is insensitive to the precise location of the matching point. Errors are 
propagated by performing the evolution on each bootstrap sample separately. They turn 
out to exceed the effect of this matching procedure by an order of magnitude. As a result 
we obtain the consistent estimates 

[ 0.1118 (10)(12)(5) fromx-T splitting 
I 0.1124 (13)(12)(15) from T' - T splitting 
We give three errors to quantify the uncertainty in the lattice scale determination: the 
statistical error, the systematic error of the NRQCD expansion and the uncertainty orig- 
inating from the sea-quark mass dependence, eq.^ permits the following conclusions: we 
have attained statistical errors on the level of the systematic uncertainties, hence they are 
not the limiting factor, even in the unquenched theory. In addition, the errors induced by 
applying nonrelativistic QCD and unphysically heavy sea quarks appear to be fairly well 
under control. 

One might be tempted to grade eq. ^as a "high precision" determination of a^: Adding 
the errors quadratically, one finds in fact an overall uncertainty of "only" ~ 1.5%, which is 
quite small compared to the errors found in recent experimental measurements of the strong 
coupling p3| . 

This conclusion, however, might be misleading. Recall that our analysis has been per- 
formed within a fixed (Wilson) discretization scheme, at a given value of the lattice cutoff 
a~^, and with an incomplete conversion prescription a^ (op ^ ), c.f. eq.^ Clearly, system- 
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atic effects which might arise from these hmitations can be taken into account properly only 
by variation of the setup. 

To estimate the size of these additional uncertainties we compare our results with those 
of ref. [^]. The latter analysis has been done at a similar value of the lattice cutoff, with 
the same conversion prescription a^^ {ap ), but within the Kogut Susskind discretization 



MS 



scheme. As their final result, the authors of ref. p quote aij^(m^) to be 0.1174(15) and 
0.1173(21) from the x — T and T' — T splitting, respectively. These values exceed our result 
by 5% or three sigmas! Thus, one concludes that the "true" systematic uncertainty is 3 to 
4 times larger than the one given in eq. ^. 

Let us discuss the possible origins of the additional uncertainty in some more detail: 
(1) Both, the errors caused by the unknown flavour dependence of C2 and the truncation 
of the perturbative series, Eq. ^ have been ignored up to now. To estimate their magnitude 
we vary the A^/-dependent part of C2 within a reasonable range, allowing values between -1 
and +1. In addition we set the coefficient of the third-order contribution to unity. The latter 
turns out to have practically no effect, whereas the variation of C2 suggests an extra 2-3% 



uncertainty. It is thus not implausible that the discrepancy in a^ is due to our ignorance 



on C2. 



(2) The difference in couplings between Wilson and staggered data is already present prior 

(2) 
to converting to the continuum renormalization scheme. We have evolved ap as it results 

from our analysis to the momentum scale used in ref. [^ and find ap (8.2 GeV) = 0.1714(25), 

a 4% difference compared to the staggered result. This deviation can be traced back directly 

to the difference in plaquette values which is larger than anticipated from the perturbative 

expansion, while the scales in both simulations are compatible within errors. This signals 

the presence of sizeable finite-a errors. 



^Of course, the same error (5%) should be added to the result of ref. |5| since, from the current stage 
of knowledge, one cannot tell which discretization scheme yields results closer to the continuum. 
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Both points constitute substantial limitations. 

To reduce the error stemming from source (1) it is of utmost importance to calculate 
perturbatively the coefficient C2{Nf) both, for Wilson and for Kogut Susskind fermions. 

A reduction of the uncertainty related to source (2) will require a much more detailed 
numerical analysis. Since the continuum limit a ^ does not exist in NRQCD, one cannot 
remove cutoff effects by extrapolation in a. Instead one has to rely on improved discretization 
schemes, which avoid sizeable cutoff effects already at finite a. The compelling test however, 
whether a given scheme really reduces cutoff effects compared to Wilson or Kogut Susskind 
discretizations can be performed only by a scaling analysis in full relativistic lattice QCD. 
Thus, in a sense, improvement of NRQCD presupposes the improvement of relativistic lattice 
QCD with respect to discretization errors. 

VI. SUMMARY AND CONCLUSION 

Heavy quarkonium bound states are potentially able to reveal the QCD coupling to high 
precision. Our main objective here has been to acquire a better understanding of the various 
sources of error when applying NRQCD techniques to carry out this program in lattice QCD. 
We have been able to reduce statistical errors to the level of systematic effects. Among the 
latter, uncertainties from the truncation of the NRQCD action and the dynamical quark 
mass dependence are found to be under good control. On the other hand, errors due to 
flavour extrapolation are more subtle to pin down but seemingly not dominant. Much 
more relevant is the choice of lattice action for the light quarks. Our analysis, using Wilson 
quarks, leads to a value of a^{mz) significantly smaller than comparable calculations based 
on staggered light quarks. This suggests that discretization errors do play an important role 
in limiting the precision of NRQCD type determinations of the strong coupling and suggests 
to base the analysis on the use of improved Wilson type actions. Needless to say, as a first 



step, one must improve the perturbative recoupling between the lattice and MS schemes. 
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TABLES 
TABLE L Simulation set-up. Where two numbers are quoted they refer to S-wave/P-wave 
correlators. 







Dynamical Wilson 




Quenched 






/3 = 5.6 , 16^ X 


32 




/? = 6.0 , 16^ X 32 


K 




0.1560 


0.1570 


0.1575 




m^^/mp 




0.833(3) 


0.758(11) 


0.686(11) 




# config 


urations 


206 


192 


203 


811/520 


^ point 


sources 


4/4 


13/8 


12/12 


4/4 



TABLE XL Fit results in lattice units. l^Si, 2'^S'i and l^Pi are obtained from two-exponential 
fits to smeared-local correlators, spin splittings are taken from ratio fits to smeared-smeared cor- 
relators. 







Nf = 2,/3 = 5.6 




Nf = 0,(3 = 6.0 


Level 


K = 0.1575 


K = 0.1570 


K = 0.1560 




l^S, 


0.3584(6) 


0.3606(6) 


0.3652(9) 


0.3438(4) 


2^ Si 


0.590(8) 


0.601(10) 


0.631(20) 


0.589(12) 


I'Pi 


0.530(5) 


0.536(6) 


0.549(7) 


0.508(6) 


l^P2 - IP 


0.0034(4) 


0.0032(4) 


0.0028(4) 


0.0032(3) 


IP - l^Pi 


0.0024(5) 


0.0021(5) 


0.0018(3) 


0.0020(4) 


IP - l3Po 


0.0099(10) 


0.0096(11) 


0.0087(10) 


0.0098(5) 


IP - l^Pi 


0.0003(4) 


0.0000(5) 


-0.0001(3) 


0.0004(3) 
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TABLE III. Results of the extrapolation in the sea-quark mass. ^P is the spin-averaged 
triplet-P state. 



splitting 


aAEo 


c 


aAE{ms/3) 


aAE{mi) 


2^50 - I'So 


0.209(21) 


1.2(7) 


0.229(10) 


0.212(19) 


2^Si - l^Si 


0.209(18) 


1.1(7) 


0.226(9) 


0.211(17) 


13P - 13^1 


0.163(9) 


0.4(3) 


0.170(5) 


0.164(8) 


2iPi - liPi 


0.152(24) 


1.8(7) 


0.181(15) 


0.156(23) 



TABLE IV. Determination of the 


lattice spacing 


from the 2^ Si - 


-I'S, 


andl^P-V 


'Si 


splittings. 




Nf = 


K = 0.1560 


K = 0.1570 




K = 0.1575 




"i-s/3 


a-i(T'-T)[GeV] 
a-i(x-T)[GeV] 


2.29(11) 
2.68(9) 


2.12(16) 

2.38(8) 


2.34(9) 
2.50(8) 




2.43(8) 
2.57(7) 




2.49(10) 

2.59(7) 



TABLE V. Comparison of quenched T and i]h energies obtained with the 0{mi,v'^) and 0{mi,v^) 
actions. Simulation parameters are (3 = 6.0, Mi, = 1.7 . 



action 


no 


l^^o 


2^50 


l35i 


2 3^1 


13^1-1 1^0 


0{mi,v^) 
0{mi,v'') 




0.3299(6) 
0.3309(4) 


0.581(12) 
0.582(12) 


0.3441(7) 
0.3438(4) 


0.587(14) 
0.589(12) 


0.01443(23) 
0.01266(8) 
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TABLE VI. Results for ap (3.41/a) extracted from the measured plaquette value. We quote 
the statistical error on the lattice scale as well as the systematic NRQCD error. 



P Nf 



M„ 



ln(iKTrt// 



fJLVl 



a 



{^/)C3.41^ 



3.41 



Mi [GeV] 



X-T 



•y^l ^ 



5.6 


2 


M,/3 


0.5570 


0.1678 




2 


Ml 


0.5546 


0.1668 


6.0 





oo 


0.5214 


0.1518 



8.84(26)(35) 8.48(33)(32) 

9.16(45)(70) 9.09(72)(71) 

9.13(34)(37) 7.82(38)(38) 



TABLE VIL Plaquette couplings at the reference scale /x = 9.0 GeV. The last two columns are 
the result of an extrapolation in the Eavour number. The first error is statistical, the second is the 
systematic error within NRQCD. 



Splitting 



aS,°^(9.0GeV) 



ag^ (9.0 GeV) 



aj?^ (9.0 GeV) 



Mg = Ms/3 



Mr, = Ml 



Mg = Ms/ 3 



Ma=Ml 



X-T 0.1525(17)(18) 0.1670(14)(19) 0.1677(24)(36) 0.1753(25)(32) 0.1764(41)(62) 
T'-T 0.1458(20)(20) 0.1650(18)(17) 0.1673(38)(37) 0.1767(33)(33) 0.1806(69)(68) 



TABLE VIIL Plaquette couplings at the reference scale /i = 9.0 GeV for each sea quark mass. 
Here, we quote only one error, covering both statistical and systematic uncertainties. 



K, 


(isRTrt/^,) 


X- 


-T 


T'- 


-T 






ag^ (9.0 GeV) 


qJ^^ (9.0 GeV) 


ag\9.0GeV) 


ag\9.0GeV) 


0.1560 


0.5698 


0.1651(26) 


0.1723(46) 


0.1598(39) 


0.1679(68) 


0.1570 


0.5716 


0.1661(25) 


0.1739(43) 


0.1631(26) 


0.1733(55) 


0.1575 


0.5725 


0.1667(22) 


0.1749(37) 


0.1642(26) 


0.1753(44) 
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TABLE IX. MS coupling at the heavy-quark thresholds and the Z mass. 
Splitting ai;5s(1.3GeV) ai^(4.1GeV) aj^(91.2GeV) 

Mq Ms/ 3 Ml Ms/3 Ml Ms/ 3 Mi 

X-T 0.316(9)(11) 0.320(14)(22) 0.2034(33) (44) 0.2050(55)(84) 0.1118(10)(12) 0.1123(16)(24) 
T'-T 0.321(12)(9) 0.335(26)(26) 0.2053(45) (37) 0.2106(95)(89) 0.1124(13)(12) 0.1139(27)(28) 
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FIGURES 
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FIG. 1. Global masses as a function of Tmin with Tmax = 30. Lines mark the selected fit value 
(solid) and its error bands (dotted). For the spin-splitting (lower right plot) smeared-local (open 
symbols) and smeared-smeared data (filled symbols) are plotted. 
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FIG. 2. Extrapolation of the splittings between the T ground state and its radial and orbital 
angular momentum excitation in the dynamical quark mass. The open triangle denotes the value 
at m. 
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FIG. 3. Bottomonium spectrum - radial and orbital angular momentum splittings. The ^Si 
ground state is constrained to match the experimental T energy. Open symbols denote quenched 
results, filled symbols Nf = 2 results at rriq = nig/S. Solid lines mark the experimental values, 
dashed lines the position of the spin-averaged ^Pj states, which turn out to be nearly identical 
with the singlet-P estimates. 
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FIG. 4. Dependence of 2^5"! — I'^Si radial splitting on the kinetic mass. 
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FIG. 5. The plaquette coupling ap as a function of the number of degenerate dynamical 
flavours. The triangles result from an extrapolation in the inverse flavour number. 
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